module PRTAllometricCarbonMod

   ! ------------------------------------------------------------------------------------
   !
   ! This module contains all of the specific functions and types for
   ! Plant Allocation and Reactive Transport Extensible Hypotheses (PARTEH)
   ! CARBON only, allometric growth hypothesis
   ! 
   ! Adapted from code originally in ED, by Rosie Fisher and Paul Moorcroft
   ! This refactor written by : Ryan Knox Apr 2018
   !
   ! ------------------------------------------------------------------------------------

  use PRTGenericMod , only  : prt_global_type
  use PRTGenericMod , only  : prt_global
  use PRTGenericMod , only  : prt_vartype
  use PRTGenericMod , only  : prt_vartypes
  use PRTGenericMod , only  : carbon12_element
  use PRTGenericMod , only  : leaf_organ
  use PRTGenericMod , only  : fnrt_organ
  use PRTGenericMod , only  : sapw_organ
  use PRTGenericMod , only  : store_organ
  use PRTGenericMod , only  : repro_organ
  use PRTGenericMod , only  : struct_organ
  use PRTGenericMod , only  : un_initialized
  use PRTGenericMod , only  : prt_carbon_allom_hyp

  use FatesAllometryMod   , only : bleaf
  use FatesAllometryMod   , only : bsap_allom
  use FatesAllometryMod   , only : bfineroot
  use FatesAllometryMod   , only : bstore_allom
  use FatesAllometryMod   , only : bdead_allom
  use FatesAllometryMod   , only : bbgw_allom
  use FatesAllometryMod   , only : bagw_allom
  use FatesAllometryMod   , only : h_allom
  use FatesAllometryMod   , only : CheckIntegratedAllometries
  use FatesAllometryMod   , only : StructureResetOfDH

  use FatesGlobals        , only : endrun => fates_endrun
  use FatesGlobals        , only : fates_log
  use shr_log_mod         , only : errMsg => shr_log_errMsg
  use FatesConstantsMod   , only : r8 => fates_r8
  use FatesConstantsMod   , only : i4 => fates_int
  use FatesIntegratorsMod , only : RKF45
  use FatesIntegratorsMod , only : Euler
  use EDPftvarcon         , only : EDPftvarcon_inst
  use FatesConstantsMod   , only : calloc_abs_error
  use FatesConstantsMod   , only : nearzero
  use FatesConstantsMod   , only : itrue
  use FatesConstantsMod   , only : years_per_day


  implicit none
  private

  ! -------------------------------------------------------------------------------------
  !
  ! Define the state variables for this specific hypothesis.
  !
  ! -------------------------------------------------------------------------------------

  integer, parameter :: leaf_c_id   = 1   ! Unique object index for leaf carbon
  integer, parameter :: fnrt_c_id   = 2   ! Unique object index for fine-root carbon
  integer, parameter :: sapw_c_id   = 3   ! Unique object index for sapwood carbon 
  integer, parameter :: store_c_id  = 4   ! Unique object index for storage carbon
  integer, parameter :: repro_c_id  = 5   ! Unique object index for reproductive carbon
  integer, parameter :: struct_c_id = 6   ! Unique object index for structural carbon
  integer, parameter :: num_vars = 6      ! THIS MUST MATCH THE LARGEST INDEX ABOVE
  
  
  ! For this hypothesis, we integrate dbh along with the other 6. Since this
  ! is a boundary condition, we do not add it to the state array, but we do want
  ! to include it with the integrator array.

  integer, parameter :: dbh_id             = 7  ! This is just used for the integrator
  integer, parameter :: n_integration_vars = 7


  ! -------------------------------------------------------------------------------------
  ! Boundary Conditions
  ! -------------------------------------------------------------------------------------
  ! Input Boundary Indices (These are public, and therefore
  !       each boundary condition across all modules must
  !       have a unique name !!!!)
  ! -------------------------------------------------------------------------------------

  integer, public, parameter :: ac_bc_inout_id_dbh   = 1   ! Plant DBH
  integer, public, parameter :: ac_bc_inout_id_netdc = 2   ! Index for the net daily C input BC
  integer, parameter         :: num_bc_inout         = 2   ! Number of in & output boundary conditions


  integer, public, parameter :: ac_bc_in_id_pft   = 1   ! Index for the PFT input BC
  integer, public, parameter :: ac_bc_in_id_ctrim = 2   ! Index for the canopy trim function
  integer, parameter         :: num_bc_in         = 2   ! Number of input boundary condition

  ! THere are no purely output boundary conditions
  integer, parameter         :: num_bc_out        = 0   ! Number of purely output boundary condtions

  ! -------------------------------------------------------------------------------------
  ! Define the size of the coorindate vector.  For this hypothesis, there is only
  ! one pool per each species x organ combination.
  ! -------------------------------------------------------------------------------------
  integer, parameter         :: icd               = 1   ! Only 1 coordinate per variable

  
  ! This is the maximum number of leaf age pools  (used for allocating scratch space)
  integer, parameter         :: max_nleafage  = 10

  ! -------------------------------------------------------------------------------------
  ! This is the core type that holds this specific
  ! plant reactive transport (PRT) module
  ! -------------------------------------------------------------------------------------


  type, public, extends(prt_vartypes) :: callom_prt_vartypes

   contains

     procedure :: DailyPRT     => DailyPRTAllometricCarbon
     procedure :: FastPRT      => FastPRTAllometricCarbon

   end type callom_prt_vartypes
   
   ! ------------------------------------------------------------------------------------
   !
   ! This next class is an extention of the base instance that maps state variables
   !      to the outside model.
   !
   ! ------------------------------------------------------------------------------------
   
   character(len=*), parameter, private :: sourcefile = __FILE__


   ! This is the instance of the mapping table and variable definitions
   ! this is only allocated once per node.  This should be read-only
   ! everywhere in the code, except for where it is populated in this init routine
   ! below.

   class(prt_global_type), public, target, allocatable :: prt_global_ac


   public :: InitPRTGlobalAllometricCarbon


contains
  
 
  subroutine InitPRTGlobalAllometricCarbon()

     ! ----------------------------------------------------------------------------------
     ! Initialize and populate the object that holds the descriptions of the variables,
     ! and contains the mappings of each variable to the pre-ordained organ
     ! and species list, and the number of boundary conditions of each 3 types.
     !
     ! This is called very early on in the call sequence of the model, and should occur
     ! before any plants start being initialized.  These mapping tables must 
     ! exist before that happens.  This initialization only happens once on each
     ! machine, and the mapping will be read-only, and a global thing. This step
     ! is not initializing the data structures bound to the plants.
     !
     ! There are two mapping tables.  One mapping table is a 2d array organized
     ! by organ and species, that contains the variable index:
     ! 
     ! prt_global%sp_organ_map
     !
     ! The other mapping table is similar, but it is a 1D array, a list of the organs.
     ! And each of these the in turn points to a list of the indices associated
     ! with that organ.  This is useful when you want to do lots of stuff to a specified
     ! organ. 
     ! 
     ! prt_global%organ_map
     !
     ! IMPORTANT NOTE:  Once this object is populated, we can use this to properly
     ! allocate the "prt_vartypes_type" objects that attached to each plant. That process
     ! is handled by generic functions, and does not need to be written in each hypothesis.
     ! 
     ! -----------------------------------------------------------------------------------

     integer :: nleafage

     allocate(prt_global_ac)
     
     ! The "state descriptor" object holds things like the names, the symbols, the units
     ! of each variable. By putting it in an object, we can loop through them when
     ! doing things like reading/writing history and restarts

     allocate(prt_global_ac%state_descriptor(num_vars))

     prt_global_ac%hyp_name = 'Allometric Carbon Only'
     
     prt_global_ac%hyp_id = prt_carbon_allom_hyp

     ! Set mapping tables to zero
     call prt_global_ac%ZeroGlobal()

     
     ! The number of leaf age classes can be determined from the parameter file,
     ! notably the size of the leaf-longevity parameter's second dimension.
     ! This is the same value in FatesInterfaceMod.F90

     nleafage = size(EDPftvarcon_inst%leaf_long,dim=2)
     
     if(nleafage>max_nleafage) then
        write(fates_log(),*) 'The allometric carbon PARTEH hypothesis'
        write(fates_log(),*) 'sets a maximum number of leaf age classes'
        write(fates_log(),*) 'used for scratch space. The model wants'
        write(fates_log(),*) 'exceed that. Simply increase max_nleafage'
        write(fates_log(),*) 'found in parteh/PRTAllometricCarbonMod.F90'
        call endrun(msg=errMsg(sourcefile, __LINE__))
     end if

     ! Register the variables. Each variable must be associated with a global identifier
     ! for an organ and species.  Leaves are a little special in that they are discretized
     ! further by age class.  Although the code works fine if this collapses to 1.

     call prt_global_ac%RegisterVarInGlobal(leaf_c_id,"Leaf Carbon","leaf_c",leaf_organ,carbon12_element,nleafage)
     call prt_global_ac%RegisterVarInGlobal(fnrt_c_id,"Fine Root Carbon","fnrt_c",fnrt_organ,carbon12_element,icd)
     call prt_global_ac%RegisterVarInGlobal(sapw_c_id,"Sapwood Carbon","sapw_c",sapw_organ,carbon12_element,icd)
     call prt_global_ac%RegisterVarInGlobal(store_c_id,"Storage Carbon","store_c",store_organ,carbon12_element,icd)
     call prt_global_ac%RegisterVarInGlobal(struct_c_id,"Structural Carbon","struct_c",struct_organ,carbon12_element,icd)
     call prt_global_ac%RegisterVarInGlobal(repro_c_id,"Reproductive Carbon","repro_c",repro_organ,carbon12_element,icd)
     
     ! Set some of the array sizes for input and output boundary conditions
     prt_global_ac%num_bc_in    = num_bc_in
     prt_global_ac%num_bc_out   = num_bc_out
     prt_global_ac%num_bc_inout = num_bc_inout
     prt_global_ac%num_vars     = num_vars

     ! Have the global generic pointer, point to this hypothesis' object
     prt_global => prt_global_ac


     return
  end subroutine InitPRTGlobalAllometricCarbon

  ! =====================================================================================
  

  subroutine DailyPRTAllometricCarbon(this)

    ! -----------------------------------------------------------------------------------
    !
    ! This is the main routine that handles allocation associated with the 1st
    ! hypothesis;  carbon only, and growth governed by allometry
    ! 
    ! This routine is explained in the technical documentation in detail.
    !
    ! Some points:
    ! 1) dbh, while not a PARTEH "state variable", is passed in from FATES (or other
    !    model), is integrated along with the mass based state variables, and then
    !    passed back to the ecosystem model. It is a "inout" style boundary condition.
    !
    ! 2) It is assumed that both growth respiration, and maintenance respiration
    !    costs have already been paid, and therefore the "carbon_balance" boundary
    !    condition is the net carbon gained by the plant over the coarse of the day.
    !    Think of "daily integrated NPP".
    ! 
    ! 3) This routine will completely spend carbon_balance if it enters as a positive 
    !    value, or replace carbon balance (using storage) if it enters as a negative value.
    !    
    ! 4) It is assumed that the ecosystem model calling this routine has ensured that
    !    the net amount of negative carbon is no greater than that which can be replaced
    !    by storage.  This routine will crash gracefully if that is not true.
    !
    ! 5) Leaves and fine-roots are given top priority, but just to replace maintenance 
    !    turnover. This can also draw from strorage.
    ! 
    ! 6) Storage is given next available carbon gain, either to push up to zero, 
    !    or to use it to top off stores.
    !
    ! 7) Third priority is then given to leaves and fine-roots again, but can only use
    !    carbon gain. Also, this transfer will attempt to get pools up to allometry.
    ! 
    ! 8) Fourth priority is to bring other live pools up to allometry, and then structure.
    ! 
    ! 9) Finally, if carbon is yet still available, it will grow all pools out concurrently
    !    including some to reproduction.
    !
    ! ----------------------------------------------------------------------------------

    
    ! The class is the only argument
    class(callom_prt_vartypes)   :: this          ! this class

    ! -----------------------------------------------------------------------------------
    ! These are local copies of the in/out boundary condition structure
    ! -----------------------------------------------------------------------------------

    real(r8),pointer :: dbh            ! Diameter at breast height [cm]
                                       ! this local will point to both in and out bc's
    real(r8),pointer :: carbon_balance ! Daily carbon balance for this cohort [kgC]

    real(r8) :: canopy_trim            ! The canopy trimming function [0-1]
    integer  :: ipft                   ! Plant Functional Type index


    real(r8) :: target_leaf_c         ! target leaf carbon [kgC]
    real(r8) :: target_fnrt_c         ! target fine-root carbon [kgC]
    real(r8) :: target_sapw_c         ! target sapwood carbon [kgC]
    real(r8) :: target_store_c        ! target storage carbon [kgC]
    real(r8) :: target_agw_c          ! target above ground carbon in woody tissues [kgC]
    real(r8) :: target_bgw_c          ! target below ground carbon in woody tissues [kgC]
    real(r8) :: target_struct_c       ! target structural carbon [kgC]

    real(r8) :: sapw_area             ! dummy var, x-section area of sapwood [m2]

    real(r8) :: leaf_below_target     ! fineroot biomass below target amount [kgC]
    real(r8) :: fnrt_below_target     ! fineroot biomass below target amount [kgC]
    real(r8) :: sapw_below_target     ! sapwood biomass below target amount [kgC]
    real(r8) :: store_below_target    ! storage biomass below target amount [kgC]
    real(r8) :: struct_below_target   ! dead (structural) biomass below target amount [kgC]
    real(r8) :: total_below_target    ! total biomass below the allometric target [kgC]

    real(r8) :: flux_adj              ! adjustment made to growth flux term to minimize error [kgC]
    real(r8) :: store_target_fraction ! ratio between storage and leaf biomass when on allometry [kgC]

    real(r8) :: leaf_c_demand         ! leaf carbon that is demanded to replace maintenance turnover [kgC]
    real(r8) :: fnrt_c_demand         ! fineroot carbon that is demanded to replace 
                                      ! maintenance turnover [kgC]
    real(r8) :: total_c_demand        ! total carbon that is demanded to replace maintenance turnover [kgC]
    logical  :: step_pass             ! Did the integration step pass?

    real(r8) :: leaf_c_flux           ! Transfer into leaves at various stages [kgC]
    real(r8) :: fnrt_c_flux           ! Transfer into fine-roots at various stages [kgC]
    real(r8) :: sapw_c_flux           ! Transfer into sapwood at various stages [kgC]
    real(r8) :: store_c_flux          ! Transfer into storage at various stages [kgC]
    real(r8) :: repro_c_flux          ! Transfer into reproduction at the final stage [kgC]
    real(r8) :: struct_c_flux         ! Transfer into structure at various stages [kgC]

    real(r8),dimension(max_nleafage) :: leaf_c0 

                                      ! Initial value of carbon used to determine net flux
    real(r8) :: fnrt_c0               ! during this routine
    real(r8) :: sapw_c0               ! ""   
    real(r8) :: store_c0              ! ""
    real(r8) :: repro_c0              ! ""
    real(r8) :: struct_c0             ! ""

    logical  :: grow_leaf             ! Are leaves at allometric target and should be grown?
    logical  :: grow_fnrt             ! Are fine-roots at allometric target and should be grown?
    logical  :: grow_sapw             ! Is sapwood at allometric target and should be grown?
    logical  :: grow_store            ! Is storage at allometric target and should be grown?

                                      ! integrator variables
    real(r8) :: deltaC                ! trial value for substep
    integer  :: ierr                  ! error flag for allometric growth step
    integer  :: nsteps                ! number of sub-steps
    integer  :: istep                 ! current substep index
    real(r8) :: totalC                ! total carbon allocated over alometric growth step
    real(r8) :: hite_out              ! dummy height variable

    integer  :: i_var                 ! index for iterating state variables
    integer  :: i_age                 ! index for iterating leaf ages
    integer  :: nleafage              ! number of leaf age classifications
    real(r8) :: leaf_age_flux         ! carbon mass flux between leaf age classification pools


    ! Integegrator variables c_pool is "mostly" carbon variables, it also includes
    ! dbh...
    ! -----------------------------------------------------------------------------------

    real(r8),dimension(n_integration_vars) :: c_pool     ! Vector of carbon pools passed to integrator
    real(r8),dimension(n_integration_vars) :: c_pool_out ! Vector of carbon pools passed back from integrator
    logical,dimension(n_integration_vars)  :: c_mask     ! Mask of active pools during integration

    integer , parameter :: max_substeps = 300            ! Maximum allowable iterations

    real(r8), parameter :: max_trunc_error = 1.0_r8      ! Maximum allowable truncation error

    integer,  parameter :: ODESolve = 2                  ! 1=RKF45,  2=Euler

    integer,  parameter :: iexp_leaf = 1                 ! index 1 is the expanding (i.e. youngest)
                                                         ! leaf age class, and therefore
                                                         ! all new allocation goes into that pool

    real(r8) ::  intgr_params(num_bc_in)                 ! The boundary conditions to this routine,
                                                         ! are pressed into an array that is also
                                                         ! passed to the integrators

    associate( & 
          leaf_c   => this%variables(leaf_c_id)%val, &
          fnrt_c   => this%variables(fnrt_c_id)%val(icd), &
          sapw_c   => this%variables(sapw_c_id)%val(icd), &
          store_c  => this%variables(store_c_id)%val(icd), &
          repro_c  => this%variables(repro_c_id)%val(icd), &
          struct_c => this%variables(struct_c_id)%val(icd))


    ! -----------------------------------------------------------------------------------
    ! 0.
    ! Copy the boundary conditions into readable local variables.
    ! We don't use pointers for bc's that ar "in" only, only "in-out" and "out"
    ! -----------------------------------------------------------------------------------

    dbh                             => this%bc_inout(ac_bc_inout_id_dbh)%rval
    carbon_balance                  => this%bc_inout(ac_bc_inout_id_netdc)%rval

    canopy_trim                     = this%bc_in(ac_bc_in_id_ctrim)%rval
    ipft                            = this%bc_in(ac_bc_in_id_pft)%ival

    intgr_params(:)                 = un_initialized
    intgr_params(ac_bc_in_id_ctrim) = this%bc_in(ac_bc_in_id_ctrim)%rval
    intgr_params(ac_bc_in_id_pft)   = real(this%bc_in(ac_bc_in_id_pft)%ival)
    
    ! -----------------------------------------------------------------------------------
    ! I. Remember the values for the state variables at the beginning of this
    ! routines. We will then use that to determine their net allocation and reactive
    ! transport flux "%net_alloc" at the end.
    ! -----------------------------------------------------------------------------------

    nleafage = prt_global%state_descriptor(leaf_c_id)%num_pos ! Number of leaf age class

    leaf_c0(1:nleafage) = leaf_c(1:nleafage)  ! Set initial leaf carbon 
    fnrt_c0 = fnrt_c                          ! Set initial fine-root carbon
    sapw_c0 = sapw_c                          ! Set initial sapwood carbon
    store_c0 = store_c                        ! Set initial storage carbon 
    repro_c0 = repro_c                        ! Set initial reproductive carbon
    struct_c0 = struct_c                      ! Set initial structural carbon

    
    ! -----------------------------------------------------------------------------------
    ! If we have more than one leaf age classification, allow
    ! some leaf biomass to transition to the older classes.  NOTE! This is not handling
    ! losses due to turnover (ie. flux from the oldest senescing class). This is only
    ! internal.
    ! (rgk 12-15-2018: Have Chonggang confirm that aging should not be restricted
    ! to evergreens)
    ! -----------------------------------------------------------------------------------

    if(nleafage>1) then
       do i_age = 1,nleafage-1
          if (EDPftvarcon_inst%leaf_long(ipft,i_age)>nearzero) then
             leaf_age_flux   = leaf_c0(i_age) * years_per_day / EDPftvarcon_inst%leaf_long(ipft,i_age)
             leaf_c(i_age)   = leaf_c(i_age) - leaf_age_flux
             leaf_c(i_age+1) = leaf_c(i_age+1) + leaf_age_flux
          end if
       end do
    end if
    

    ! -----------------------------------------------------------------------------------
    ! II. Calculate target size of the biomass compartment for a given dbh.   
    ! -----------------------------------------------------------------------------------
    
    ! Target sapwood biomass and deriv. according to allometry and trimming [kgC, kgC/cm]
    call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c)

    ! Target total above ground deriv. biomass in woody/fibrous tissues  [kgC, kgC/cm]
    call bagw_allom(dbh,ipft,target_agw_c)

    ! Target total below ground deriv. biomass in woody/fibrous tissues [kgC, kgC/cm] 
    call bbgw_allom(dbh,ipft,target_bgw_c)

    ! Target total dead (structrual) biomass and deriv. [kgC, kgC/cm]
    call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c)


    ! ------------------------------------------------------------------------------------
    ! If structure is larger than target, then we need to correct some integration errors
    ! by slightly increasing dbh to match it.
    ! For grasses, if leaf biomass is larger than target, then we reset dbh to match
    ! -----------------------------------------------------------------------------------
    if( (( struct_c - target_struct_c ) > calloc_abs_error) .and. &
          (EDPftvarcon_inst%woody(ipft) == itrue) ) then

       call StructureResetOfDH( struct_c, ipft, &
             canopy_trim, dbh, hite_out )

       ! Target sapwood biomass and deriv. according to allometry and trimming [kgC, kgC/cm]
       call bsap_allom(dbh,ipft,canopy_trim,sapw_area,target_sapw_c)

       ! Target total above ground deriv. biomass in woody/fibrous tissues  [kgC, kgC/cm]
       call bagw_allom(dbh,ipft,target_agw_c)
       
       ! Target total below ground deriv. biomass in woody/fibrous tissues [kgC, kgC/cm] 
       call bbgw_allom(dbh,ipft,target_bgw_c)
       
       ! Target total dead (structrual) biomass and deriv. [kgC, kgC/cm]
       call bdead_allom( target_agw_c, target_bgw_c, target_sapw_c, ipft, target_struct_c)

    end if
    
    ! Target leaf biomass according to allometry and trimming
    call bleaf(dbh,ipft,canopy_trim,target_leaf_c)

    ! Target fine-root biomass and deriv. according to allometry and trimming [kgC, kgC/cm]
    call bfineroot(dbh,ipft,canopy_trim,target_fnrt_c)

    ! Target storage carbon [kgC,kgC/cm]
    call bstore_allom(dbh,ipft,canopy_trim,target_store_c)

    


    ! -----------------------------------------------------------------------------------
    ! III.  Prioritize some amount of carbon to replace leaf/root turnover
    !         Make sure it isnt a negative payment, and either pay what is available
    !         or forcefully pay from storage. 
    ! -----------------------------------------------------------------------------------
    
    if( EDPftvarcon_inst%evergreen(ipft) ==1 ) then
       leaf_c_demand   = max(0.0_r8, &
             EDPftvarcon_inst%leaf_stor_priority(ipft)*sum(this%variables(leaf_c_id)%turnover(:)))
    else
       leaf_c_demand   = 0.0_r8
    end if
    
    fnrt_c_demand   = max(0.0_r8, &
          EDPftvarcon_inst%leaf_stor_priority(ipft)*this%variables(fnrt_c_id)%turnover(icd))

    total_c_demand = leaf_c_demand + fnrt_c_demand
    
    if (total_c_demand> nearzero ) then

       ! We pay this even if we don't have the carbon
       ! Just don't pay so much carbon that storage+carbon_balance can't pay for it

       leaf_c_flux = min(leaf_c_demand, &
                         max(0.0_r8,(store_c+carbon_balance)* &
                         (leaf_c_demand/total_c_demand)))
       
       ! Add carbon to the youngest age pool (i.e iexp_leaf = index 1)
       carbon_balance    = carbon_balance   - leaf_c_flux
       leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux

       ! If we are testing b4b, then we pay this even if we don't have the carbon
       fnrt_c_flux = min(fnrt_c_demand, &
                         max(0.0_r8, (store_c+carbon_balance)* &
                         (fnrt_c_demand/total_c_demand)))

       carbon_balance = carbon_balance - fnrt_c_flux
       fnrt_c         = fnrt_c + fnrt_c_flux

    end if

    ! -----------------------------------------------------------------------------------
    ! IV. if carbon balance is negative, re-coup the losses from storage
    !       if it is positive, give some love to storage carbon
    ! -----------------------------------------------------------------------------------

    if( carbon_balance < 0.0_r8 ) then
       
       store_c_flux           = carbon_balance
       carbon_balance         = carbon_balance - store_c_flux
       store_c                = store_c + store_c_flux

    else

       store_below_target     = max(target_store_c - store_c,0.0_r8)
       store_target_fraction  = max(0.0_r8, store_c/target_store_c )

       store_c_flux           = min(store_below_target,carbon_balance * &
                                max(exp(-1.*store_target_fraction**4._r8) - exp( -1.0_r8 ),0.0_r8))

       carbon_balance         = carbon_balance - store_c_flux
       store_c                = store_c + store_c_flux

    end if

    ! -----------------------------------------------------------------------------------
    ! V.  If carbon is still available, prioritize some allocation to replace
    !        the rest of the leaf/fineroot deficit
    !        carbon balance is guaranteed to be >=0 beyond this point
    ! -----------------------------------------------------------------------------------
    
    leaf_c_demand   = max(0.0_r8,(target_leaf_c - sum(leaf_c(1:nleafage))))
    fnrt_c_demand   = max(0.0_r8,(target_fnrt_c - fnrt_c))

    total_c_demand = leaf_c_demand + fnrt_c_demand
    
    if( (carbon_balance > nearzero ) .and. (total_c_demand>nearzero)) then

       leaf_c_flux    = min(leaf_c_demand, &
                        carbon_balance*(leaf_c_demand/total_c_demand))
       carbon_balance = carbon_balance - leaf_c_flux
       leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux
       
       fnrt_c_flux    = min(fnrt_c_demand, &
                            carbon_balance*(fnrt_c_demand/total_c_demand))
       carbon_balance = carbon_balance - fnrt_c_flux
       fnrt_c         = fnrt_c + fnrt_c_flux

    end if

    ! -----------------------------------------------------------------------------------
    ! VI.  If carbon is still available, we try to push all live 
    !        pools back towards allometry. But only upwards, if fusion happened
    !        to generate some pools above allometric target, don't reduce the pool,
    !        just ignore it until the rest of the plant grows to meet it.
    ! -----------------------------------------------------------------------------------
    if( carbon_balance > nearzero ) then

       leaf_below_target  = max(target_leaf_c - sum(leaf_c(1:nleafage)),0.0_r8)
       fnrt_below_target  = max(target_fnrt_c - fnrt_c,0.0_r8)
       sapw_below_target  = max(target_sapw_c - sapw_c,0.0_r8)
       store_below_target = max(target_store_c - store_c,0.0_r8)
       
       total_below_target = leaf_below_target + fnrt_below_target + &
                            sapw_below_target + store_below_target
    
       if ( total_below_target > nearzero ) then
          
          if( total_below_target > carbon_balance) then
             leaf_c_flux  = carbon_balance * leaf_below_target/total_below_target
             fnrt_c_flux = carbon_balance * fnrt_below_target/total_below_target
             sapw_c_flux  = carbon_balance * sapw_below_target/total_below_target
             store_c_flux = carbon_balance * store_below_target/total_below_target
          else
             leaf_c_flux  = leaf_below_target
             fnrt_c_flux = fnrt_below_target
             sapw_c_flux = sapw_below_target
             store_c_flux = store_below_target
          end if

          carbon_balance               = carbon_balance - leaf_c_flux
          leaf_c(iexp_leaf)            = leaf_c(iexp_leaf) + leaf_c_flux
          
          carbon_balance               = carbon_balance - fnrt_c_flux
          fnrt_c                       = fnrt_c + fnrt_c_flux
          
          carbon_balance               = carbon_balance - sapw_c_flux
          sapw_c                       = sapw_c + sapw_c_flux
          
          carbon_balance               = carbon_balance - store_c_flux
          store_c                      = store_c  +  store_c_flux
          
       end if
    end if
    
    ! -----------------------------------------------------------------------------------
    ! VII.  If carbon is still available, replenish the structural pool to get
    !           back on allometry
    ! -----------------------------------------------------------------------------------

    if( carbon_balance > nearzero ) then
    
       struct_below_target  = max(target_struct_c - struct_c ,0.0_r8)
    
       if ( struct_below_target > 0.0_r8) then
       
          struct_c_flux          = min(carbon_balance,struct_below_target)
          carbon_balance         = carbon_balance - struct_c_flux
          struct_c               = struct_c + struct_c_flux
          
       end if

    end if
    
    ! -----------------------------------------------------------------------------------
    ! VIII.  If carbon is yet still available ...
    !        Our pools are now either on allometry or above (from fusion).
    !        We we can increment those pools at or below,
    !        including structure and reproduction according to their rates
    !        Use an adaptive euler integration. If the error is not nominal,
    !        the carbon balance sub-step (deltaC) will be halved and tried again
    !
    ! Note that we compare against calloc_abs_error here because it is possible
    ! that all the carbon was effectively used up, but a miniscule amount
    ! remains due to numerical precision (ie -20 or so), so even though
    ! the plant has not been brought to be "on allometry", it thinks it has carbon
    ! left to allocate, and thus it must be on allometry when its not.
    ! -----------------------------------------------------------------------------------
    
    if( carbon_balance > calloc_abs_error ) then
       
       ! This routine checks that actual carbon is not below that targets. It does
       ! allow actual pools to be above the target, and in these cases, it sends
       ! a false on the "grow_<>" flag, allowing the plant to grow into these pools.
       ! It also checks to make sure that structural biomass is not above the target.
       if ( EDPftvarcon_inst%woody(ipft) == itrue ) then


          if( (target_store_c - store_c)>calloc_abs_error) then
             write(fates_log(),*) 'storage is not on-allometry at the growth step'
             write(fates_log(),*) 'exiting'
             write(fates_log(),*) 'cbal: ',carbon_balance
             write(fates_log(),*) 'near-zero',nearzero
             write(fates_log(),*) 'store_c: ',store_c
             write(fates_log(),*) 'target c: ',target_store_c
             write(fates_log(),*) 'store_c0:', store_c0
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if


          call TargetAllometryCheck(sum(leaf_c(1:nleafage)), fnrt_c, sapw_c, &
                                    store_c, struct_c,       &
                                    target_leaf_c, target_fnrt_c, &
                                    target_sapw_c, target_store_c, target_struct_c, &
                                    grow_leaf, grow_fnrt, grow_sapw, grow_store)
       else ! for grasses 
          grow_leaf  = .true.
          grow_fnrt = .true.
          grow_sapw  = .true.
          grow_store = .true.
       end if

       ! --------------------------------------------------------------------------------
       ! The numerical integration of growth requires that the instantaneous state
       ! variables are passed in as an array.  We call it "c_pool".
       !
       ! Initialize the adaptive integrator arrays and flags
       ! --------------------------------------------------------------------------------

       ierr             = 1
       totalC           = carbon_balance
       nsteps           = 0
       
       c_pool(:) = 0.0_r8                        ! Zero state variable array
       c_mask(:) = .false.                       ! This mask tells the integrator
                                                 ! which indices are active. Its possible
                                                 ! that due to fusion, or previous numerical
                                                 ! truncation errors, that one of these pools
                                                 ! may be larger than its target! We check
                                                 ! this, and if true, then we flag that
                                                 ! pool to be ignored. c_mask(i) = .false.
                                                 ! For grasses, since they don't grow very 
                                                 ! large and thus won't accumulate such large
                                                 ! errors, we always mask as true.

       c_pool(leaf_c_id)   = sum(leaf_c(1:nleafage))
       c_pool(fnrt_c_id)   = fnrt_c
       c_pool(sapw_c_id)   = sapw_c
       c_pool(store_c_id)  = store_c
       c_pool(struct_c_id) = struct_c
       c_pool(repro_c_id)  = repro_c
       c_pool(dbh_id)      = dbh
       
       c_mask(leaf_c_id)   = grow_leaf
       c_mask(fnrt_c_id)   = grow_fnrt
       c_mask(sapw_c_id)   = grow_sapw
       c_mask(store_c_id)  = grow_store
       c_mask(struct_c_id) = .true.                ! Always increment dead on growth step
       c_mask(repro_c_id)  = .true.                ! Always calculate reproduction on growth
       c_mask(dbh_id)      = .true.                ! Always increment dbh on growth step
       

       ! When using the Euler method, we keep things simple.  We always try
       ! to make the first integration step to span the entirety of the integration
       ! window for the independent variable (available carbon)

       if(ODESolve == 2) then
          this%ode_opt_step = totalC
       end if
       
       do while( ierr .ne. 0 )
          
          deltaC = min(totalC,this%ode_opt_step)
          if(ODESolve == 1) then
             call RKF45(AllomCGrowthDeriv,c_pool,c_mask,deltaC,totalC, &
                   max_trunc_error,intgr_params,c_pool_out,this%ode_opt_step,step_pass)
             
          elseif(ODESolve == 2) then
             call Euler(AllomCGrowthDeriv,c_pool,c_mask,deltaC,totalC,intgr_params,c_pool_out)
             !  step_pass = .true.
             
             ! When integrating along the allometric curve, we have the luxury of perfect
             ! hindsite.  Ie, after we have made our step, we can see if the amount
             ! of each carbon we have matches the target associated with the new dbh.
             ! The following call evaluates how close we are to the allometically defined
             ! targets. If we are too far (governed by max_trunc_error), then we
             ! pass back the pass/fail flag (step_pass) as false.  If false, then
             ! we halve the step-size, and then retry.  If that step was fine, then
             ! we remember the current step size as a good next guess.
             
             call CheckIntegratedAllometries(c_pool_out(dbh_id),ipft,canopy_trim,  &
                   c_pool_out(leaf_c_id), c_pool_out(fnrt_c_id), c_pool_out(sapw_c_id), &
                   c_pool_out(store_c_id), c_pool_out(struct_c_id), &
                   c_mask(leaf_c_id), c_mask(fnrt_c_id), c_mask(sapw_c_id), &
                   c_mask(store_c_id),c_mask(struct_c_id),  max_trunc_error, step_pass)
             if(step_pass)  then
                this%ode_opt_step = deltaC
             else
                this%ode_opt_step = 0.5*deltaC
             end if
          else
             write(fates_log(),*) 'An integrator was chosen that does not exist'
             write(fates_log(),*) 'ODESolve = ',ODESolve
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if

          nsteps = nsteps + 1
          
          if (step_pass) then ! If true, then step is accepted
             totalC    = totalC - deltaC
             c_pool(:) = c_pool_out(:)
          end if
          
          if(nsteps > max_substeps ) then
             write(fates_log(),*) 'Plant Growth Integrator could not find'
             write(fates_log(),*) 'a solution in less than ',max_substeps,' tries'
             write(fates_log(),*) 'Aborting'
             write(fates_log(),*) 'carbon_balance',carbon_balance
             write(fates_log(),*) 'deltaC',deltaC
             write(fates_log(),*) 'totalC',totalC
             write(fates_log(),*) 'leaf:',grow_leaf,target_leaf_c,target_leaf_c - sum(leaf_c(:))
             write(fates_log(),*) 'fnrt:',grow_fnrt,target_fnrt_c,target_fnrt_c - fnrt_c
             write(fates_log(),*) 'sap:',grow_sapw,target_sapw_c, target_sapw_c - sapw_c
             write(fates_log(),*) 'store:',grow_store,target_store_c,target_store_c - store_c
             write(fates_log(),*) 'dead:',target_struct_c,target_struct_c - struct_c
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if

          !
          ! TotalC should eventually be whittled down to near zero.
          ! The solvers are not perfect, so we can't expect it to be perfectly zero.
          ! Note that calloc_abs_error is 1e-9, which is really small (1 microgram of carbon)
          ! yet also six orders of magnitude greater than typical rounding errors (~1e-15).
  
          ! At that point, update the actual states
          ! --------------------------------------------------------------------------------
          if( (totalC < calloc_abs_error) .and. (step_pass) )then

             ierr           = 0
             leaf_c_flux    = c_pool(leaf_c_id)   - sum(leaf_c(1:nleafage))
             fnrt_c_flux    = c_pool(fnrt_c_id)   - fnrt_c
             sapw_c_flux    = c_pool(sapw_c_id)   - sapw_c
             store_c_flux   = c_pool(store_c_id)  - store_c
             struct_c_flux  = c_pool(struct_c_id) - struct_c
             repro_c_flux   = c_pool(repro_c_id)  - repro_c
             
             ! Make an adjustment to flux partitions to make it match remaining c balance
             flux_adj       = carbon_balance/(leaf_c_flux+fnrt_c_flux+sapw_c_flux + &
                                              store_c_flux+struct_c_flux+repro_c_flux)

             
             leaf_c_flux    = leaf_c_flux*flux_adj
             fnrt_c_flux    = fnrt_c_flux*flux_adj
             sapw_c_flux    = sapw_c_flux*flux_adj
             store_c_flux   = store_c_flux*flux_adj
             struct_c_flux  = struct_c_flux*flux_adj
             repro_c_flux   = repro_c_flux*flux_adj
             
             carbon_balance    = carbon_balance - leaf_c_flux
             leaf_c(iexp_leaf) = leaf_c(iexp_leaf) + leaf_c_flux
             
             carbon_balance = carbon_balance - fnrt_c_flux
             fnrt_c         = fnrt_c + fnrt_c_flux
             
             carbon_balance = carbon_balance - sapw_c_flux
             sapw_c         = sapw_c + sapw_c_flux
             
             carbon_balance = carbon_balance - store_c_flux
             store_c        = store_c + store_c_flux
             
             carbon_balance = carbon_balance - struct_c_flux
             struct_c       = struct_c + struct_c_flux
             
             carbon_balance = carbon_balance - repro_c_flux
             repro_c        = repro_c  + repro_c_flux
             
             dbh            = c_pool(dbh_id)

             ! THESE HAVE TO BE SET OUTSIDE OF THIS ROUTINE
             !!          cohort%seed_prod = cohort%seed_prod + brepro_flux / hlm_freq_day
             !!          cohort%dhdt      = (h_sub-cohort%hite)/hlm_freq_day
             !!          cohort%ddbhdt    = (dbh_sub-dbh_in)/hlm_freq_day
             
             if( abs(carbon_balance)>calloc_abs_error ) then
                write(fates_log(),*) 'carbon conservation error while integrating pools'
                write(fates_log(),*) 'along alometric curve'
                write(fates_log(),*) 'carbon_balance = ',carbon_balance,totalC
                write(fates_log(),*) 'exiting'
                call endrun(msg=errMsg(sourcefile, __LINE__))
             end if
             
          end if
       end do
    end if

    ! Track the net allocations and transport from this routine

    do i_age = 1,nleafage
       this%variables(leaf_c_id)%net_alloc(i_age) = &
             this%variables(leaf_c_id)%net_alloc(i_age) + &
             (leaf_c(i_age) - leaf_c0(i_age))
    end do

    this%variables(fnrt_c_id)%net_alloc(icd) = &
         this%variables(fnrt_c_id)%net_alloc(icd) + (fnrt_c - fnrt_c0)
    
    this%variables(sapw_c_id)%net_alloc(icd) = &
         this%variables(sapw_c_id)%net_alloc(icd) + (sapw_c - sapw_c0)
    
    this%variables(store_c_id)%net_alloc(icd) = &
         this%variables(store_c_id)%net_alloc(icd) + (store_c - store_c0)
    
    this%variables(repro_c_id)%net_alloc(icd) = &
         this%variables(repro_c_id)%net_alloc(icd) + (repro_c - repro_c0)
    
    this%variables(struct_c_id)%net_alloc(icd) = &
         this%variables(struct_c_id)%net_alloc(icd) + (struct_c - struct_c0)



   end associate
  
   return
  end subroutine DailyPRTAllometricCarbon
  
  ! =====================================================================================
  
  function AllomCGrowthDeriv(c_pools,c_mask,cbalance,intgr_params) result(dCdx)

      ! ---------------------------------------------------------------------------------
      ! This function calculates the derivatives for the carbon pools
      ! relative to the amount of carbon balance.  This function is based completely
      ! off of allometry, and assumes that there are no other species (ie nutrients) that
      ! govern allocation.
      ! ---------------------------------------------------------------------------------
      
      ! Arguments
      real(r8),intent(in), dimension(:) :: c_pools      ! Vector of carbon pools
                                                        ! dbh,leaf,root,sap,store,dead
      logical,intent(in), dimension(:)  :: c_mask       ! logical mask of active pools
                                                        ! some may be turned off
      real(r8),intent(in)               :: cbalance     ! The carbon balance of the
                                                        ! partial step (independant var)
                                             
      real(r8), intent(in),dimension(:) :: intgr_params  ! Generic Array used to pass 
                                                         ! parameters into this function


      ! Return Value 
      ! Change in carbon (each pool) per change in total allocatable carbon (kgC/kgC)
      real(r8),dimension(lbound(c_pools,dim=1):ubound(c_pools,dim=1)) :: dCdx 

      ! locals
      integer  :: ipft       ! PFT index
      real(r8) :: canopy_trim    ! Canopy trimming function (boundary condition [0-1]
      real(r8) :: ct_leaf    ! target leaf biomass, dummy var (kgC)
      real(r8) :: ct_fnrt   ! target fine-root biomass, dummy var (kgC)
      real(r8) :: ct_sap     ! target sapwood biomass, dummy var (kgC)
      real(r8) :: ct_agw     ! target aboveground wood, dummy var (kgC)
      real(r8) :: ct_bgw     ! target belowground wood, dummy var (kgC)
      real(r8) :: ct_store   ! target storage, dummy var (kgC)
      real(r8) :: ct_dead    ! target structural biomas, dummy var (kgC)
      real(r8) :: sapw_area      ! dummy sapwood area
      real(r8) :: ct_dleafdd     ! target leaf biomass derivative wrt diameter, (kgC/cm)
      real(r8) :: ct_dfnrtdd     ! target fine-root biomass derivative wrt diameter, (kgC/cm)
      real(r8) :: ct_dsapdd      ! target sapwood biomass derivative wrt diameter, (kgC/cm)
      real(r8) :: ct_dagwdd      ! target AG wood biomass derivative wrt diameter, (kgC/cm)
      real(r8) :: ct_dbgwdd      ! target BG wood biomass derivative wrt diameter, (kgC/cm)
      real(r8) :: ct_dstoredd    ! target storage biomass derivative wrt diameter, (kgC/cm)
      real(r8) :: ct_ddeaddd     ! target structural biomass derivative wrt diameter, (kgC/cm)
      real(r8) :: ct_dtotaldd    ! target total (not reproductive) biomass derivative wrt diameter, (kgC/cm)
      real(r8) :: repro_fraction ! fraction of carbon balance directed towards reproduction (kgC/kgC)


      associate( dbh    => c_pools(dbh_id), &
                 cleaf  => c_pools(leaf_c_id), &
                 cfnrt  => c_pools(fnrt_c_id), &
                 csap   => c_pools(sapw_c_id), &
                 cstore => c_pools(store_c_id), &
                 cdead  => c_pools(struct_c_id), &
                 crepro => c_pools(repro_c_id), &    ! Unused (memoryless)
                 mask_dbh  => c_mask(dbh_id), &    ! Unused (dbh always grows)
                 mask_leaf => c_mask(leaf_c_id), &  
                 mask_fnrt => c_mask(fnrt_c_id), &
                 mask_sap  => c_mask(sapw_c_id), &
                 mask_store => c_mask(store_c_id), &
                 mask_dead  => c_mask(struct_c_id),  & ! Unused (dead always grows)
                 mask_repro => c_mask(repro_c_id) )

        canopy_trim = intgr_params(ac_bc_in_id_ctrim)
        ipft        = int(intgr_params(ac_bc_in_id_pft))

        call bleaf(dbh,ipft,canopy_trim,ct_leaf,ct_dleafdd)
        call bfineroot(dbh,ipft,canopy_trim,ct_fnrt,ct_dfnrtdd)
        call bsap_allom(dbh,ipft,canopy_trim,sapw_area,ct_sap,ct_dsapdd)

        call bagw_allom(dbh,ipft,ct_agw,ct_dagwdd)
        call bbgw_allom(dbh,ipft,ct_bgw,ct_dbgwdd)
        call bdead_allom(ct_agw,ct_bgw, ct_sap, ipft, ct_dead, &
                         ct_dagwdd, ct_dbgwdd, ct_dsapdd, ct_ddeaddd)
        call bstore_allom(dbh,ipft,canopy_trim,ct_store,ct_dstoredd)
        
        ! fraction of carbon going towards reproduction
        if (dbh <= EDPftvarcon_inst%dbh_repro_threshold(ipft)) then ! cap on leaf biomass
           repro_fraction = EDPftvarcon_inst%seed_alloc(ipft)
        else
           repro_fraction = EDPftvarcon_inst%seed_alloc(ipft) + EDPftvarcon_inst%seed_alloc_mature(ipft)
        end if

        dCdx = 0.0_r8

        ct_dtotaldd = ct_ddeaddd
        if (mask_leaf)  ct_dtotaldd = ct_dtotaldd + ct_dleafdd 
        if (mask_fnrt) ct_dtotaldd = ct_dtotaldd + ct_dfnrtdd 
        if (mask_sap)   ct_dtotaldd = ct_dtotaldd + ct_dsapdd
        if (mask_store) ct_dtotaldd = ct_dtotaldd + ct_dstoredd

        ! It is possible that with some asymptotic, or hard
        ! capped allometries, that all growth rates reach zero.
        ! In this case, if there is carbon, give it to reproduction

        if(ct_dtotaldd<=tiny(ct_dtotaldd))then

           dCdx(struct_c_id)  = 0.0_r8
           dCdx(dbh_id)    = 0.0_r8
           dCdx(leaf_c_id)  = 0.0_r8
           dCdx(fnrt_c_id) = 0.0_r8
           dCdx(sapw_c_id)   = 0.0_r8
           dCdx(store_c_id) = 0.0_r8
           dCdx(repro_c_id) = 1.0_r8

        else

           dCdx(struct_c_id) = (ct_ddeaddd/ct_dtotaldd)*(1.0_r8-repro_fraction)   
           dCdx(dbh_id)   = (1.0_r8/ct_dtotaldd)*(1.0_r8-repro_fraction)
        
           if (mask_leaf) then
              dCdx(leaf_c_id) = (ct_dleafdd/ct_dtotaldd)*(1.0_r8-repro_fraction)
           else
              dCdx(leaf_c_id) = 0.0_r8
           end if
           
           if (mask_fnrt) then
              dCdx(fnrt_c_id) = (ct_dfnrtdd/ct_dtotaldd)*(1.0_r8-repro_fraction)
           else
              dCdx(fnrt_c_id) = 0.0_r8
           end if
           
           if (mask_sap) then
              dCdx(sapw_c_id) = (ct_dsapdd/ct_dtotaldd)*(1.0_r8-repro_fraction)
           else
              dCdx(sapw_c_id) = 0.0_r8
           end if
           
           if (mask_store) then
              dCdx(store_c_id) = (ct_dstoredd/ct_dtotaldd)*(1.0_r8-repro_fraction)
           else
              dCdx(store_c_id) = 0.0_r8
           end if
           
           dCdx(repro_c_id) = repro_fraction

        end if
        
      end associate

      return
   end function AllomCGrowthDeriv

   ! ====================================================================================

   subroutine TargetAllometryCheck(bleaf,bfroot,bsap,bstore,bdead, &
                                   bt_leaf,bt_froot,bt_sap,bt_store,bt_dead, &
                                   grow_leaf,grow_froot,grow_sapw,grow_store)

     ! Arguments
     real(r8),intent(in) :: bleaf   !actual
     real(r8),intent(in) :: bfroot
     real(r8),intent(in) :: bsap
     real(r8),intent(in) :: bstore
     real(r8),intent(in) :: bdead
     real(r8),intent(in) :: bt_leaf   !target
     real(r8),intent(in) :: bt_froot
     real(r8),intent(in) :: bt_sap
     real(r8),intent(in) :: bt_store
     real(r8),intent(in) :: bt_dead
     logical,intent(out) :: grow_leaf  !growth flag
     logical,intent(out) :: grow_froot
     logical,intent(out) :: grow_sapw
     logical,intent(out) :: grow_store
       
     if( (bt_leaf - bleaf)>calloc_abs_error) then
        write(fates_log(),*) 'leaves are not on-allometry at the growth step'
        write(fates_log(),*) 'exiting',bleaf,bt_leaf
        call endrun(msg=errMsg(sourcefile, __LINE__))
     elseif( (bleaf - bt_leaf)>calloc_abs_error) then
        ! leaf is above allometry, ignore
        grow_leaf = .false.
     else
        grow_leaf = .true.
     end if
     
     if( (bt_froot - bfroot)>calloc_abs_error) then
        write(fates_log(),*) 'fineroots are not on-allometry at the growth step'
        write(fates_log(),*) 'exiting',bfroot, bt_froot
        call endrun(msg=errMsg(sourcefile, __LINE__))
     elseif( ( bfroot-bt_froot)>calloc_abs_error ) then
        grow_froot = .false.
     else
        grow_froot = .true.
     end if
     
     if( (bt_sap - bsap)>calloc_abs_error) then
        write(fates_log(),*) 'sapwood is not on-allometry at the growth step'
        write(fates_log(),*) 'exiting',bsap, bt_sap
        call endrun(msg=errMsg(sourcefile, __LINE__))
     elseif( ( bsap-bt_sap)>calloc_abs_error ) then
        grow_sapw = .false.
     else
        grow_sapw = .true.
     end if
     
     if( (bt_store - bstore)>calloc_abs_error) then
        write(fates_log(),*) 'storage is not on-allometry at the growth step'
        write(fates_log(),*) 'exiting',bstore,bt_store
        call endrun(msg=errMsg(sourcefile, __LINE__))
     elseif( ( bstore-bt_store)>calloc_abs_error ) then
        grow_store = .false.
     else
        grow_store = .true.
     end if
     
     if( (bt_dead - bdead)>calloc_abs_error) then
        write(fates_log(),*) 'structure not on-allometry at the growth step'
        write(fates_log(),*) 'exiting',bdead,bt_dead
        call endrun(msg=errMsg(sourcefile, __LINE__))
     end if
   end subroutine TargetAllometryCheck

   ! =====================================================================================

   subroutine FastPRTAllometricCarbon(this)
      
      implicit none
      class(callom_prt_vartypes) :: this     ! this class
      
      ! This routine does nothing, because in the carbon only allometric RT model
      ! we currently don't have any fast-timestep processes
      ! Think of this as a stub.


      return
   end subroutine FastPRTAllometricCarbon


end module PRTAllometricCarbonMod
  
